Mathematics of Statistical Parallax and the Local Distance Scale 
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ABSTRACT 

We present a mathematical analysis of the statistical parallax method. The 
method yields physical insight into the maximum-likelihood determinations of 
the luminosity and velocity distribution and enables us to conduct a vigorous 
Monte Carlo investigation into various systematic effects. We apply our analytic 
formalism to the RR Lyrae sample of Layden et al. The velocity distribution 
of RR Lyrae stars is highly non-Gaussian, with kurtoses K n = 2.04, K e = 3.22 
and K z = 4.28 in the three principal directions, but this has almost no effect on 
either the best fit or the uncertainty of the luminosity determination. Indeed, 
our principal result is that the statistical parallax method is extremely robust 
in the face of all systematic effects that we considered. 

Our analysis, applied to the Layden et al. RR Lyrae sample, strictly confirms 
the majority of their results. The mean RR Lyrae absolute magnitude is 
M v = 0.75 ± 0.13 at the mean metallicity of the sample ([Fe/H]) = -1.61, 
compared to My = 0.71 ± 0.12 obtained by Layden et al. Most of the difference 
is due to Malmquist bias which was not considered in previous studies. We 
also analyze a semi-independent non-kinematically selected sample of stars 
with metallicities at the [Fe/H]< —1.5 taken from Layden et al. and Beers & 
Sommer-Larsen and obtain My = 0.79 ±0.12 at ([Fe/H]) = -1.79. Additionally, 
this analysis yields measurements of the radial bulk motion (4 ± lOkms -1 ) and 
vertical bulk motion (0 ± 6kms _1 ) of the halo relative to the Local Standard of 
Rest. 

Subject headings: distance scale — Galaxy: kinematics and dynamics- 
methods: analytical, statistical — stars: variables: RR Lyrae 
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1. Introduction 

The discrepancy between the Cepheid and RR Lyrae distance estimates to the Large 
Magellanic Cloud (LMC) is a long-standing problem (e.g., van den Bergh 1995). There are 
two possible sources of disagreement between the Cepheid and RR Lyrae distance scales: 
either the compared populations of Cepheids or RR Lyrae stars in the Milky Way and LMC 
are different or there are substantial errors in the Cepheid or RR Lyrae calibration. In this 
paper we investigate the statistical and systematic errors associated with the statistical 
parallax method, which is one of the main calibration methods for RR Lyrae stars. We 
prove that the statistical parallax method is very robust. Consequently, miscalibration of 
the absolute magnitude of RR Lyrae stars is not likely to be the solution to the distance 
scale problem and it is fair to state that RR Lyrae stars constitute reliable distance 
indicators. We conclude that either the Cepheid distance scale is incorrect or the currently 
available data do not allow for a reliable comparison between RR Lyrae stars in the Milky 
Way and LMC. 

Historically, the methods of secular parallax and classical statistical parallax (e.g., 
Trumpler & Weaver 1962) have been used to expand the method of trigonometric parallax 
to more distant objects by increasing the baseline of the measurement. Secular parallax 
takes advantage of the fact that the Sun moves at a speed W ~ 20kms~ 1 relative to disk 
populations and W ~ 200 km s -1 relative to halo populations, corresponding to 4 AU per 
year and 40 AU per year, respectively. After a time t of several decades, the accumulated 
baseline Wt is several orders of magnitude larger than the Earth-Sun baseline used in 
trigonometric parallax. If the speed of the individual star (rather than of the population 
as a whole) were known to be W, then one could measure the distance with a precision 
better by a factor Wt /AU than using trigonometric parallax for the same astrometric 
accuracy. The drawback is that the stars have a dispersion a, so that even in the limit of 
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perfect measurements, the precision of the distance determination is inversely proportional 
to the "Mach number", k = W/cr, which is typically of order 1. By applying this method 
to a sample of N stars, one can beat down the noise and achieve a precision oc iV -1 / 2 /? -1 . 
The velocity dispersion itself can be used to obtain an independent measurement of the 
distance using the method of classical statistical parallax. That is, the distance scale of 
a class of stars can be estimated by forcing their radial velocities and proper motions to 
reproduce the same velocity dispersion. The precision of this method is oc iV -1 / 2 . The two 
determinations can be combined to achieve a joint precision oc [N(l + g^c 2 )]" 1 / 2 , where 
g ~ 1/6 is a geometrical factor that we derive below. Nevertheless, because the method is 
limited by the number of stars rather than the precision of the measurements, its usefulness 
has been restricted to ever more distant classes of stars as astrometry has improved. 

Only recently (e.g., Murray 1983; Hawley et al. 1986; Strugnell, Reid, & Murray 1986) 
was it fully realized that secular parallax and classical statistical parallax are actually two 
aspects of the same generalized method, now simply called "statistical parallax" . Secular 
parallax is based on forcing equality between the three first moments of the velocity 
distribution (the bulk motion W) as determined from radial velocity and proper motion 
measurements, while classical statistical parallax is based on forcing equality of the six 
second moments (the six independent components of the velocity covariance matrix C^). In 
the modern combined version of statistical parallax one simply determines ten parameters 
simultaneously by applying maximum likelihood. The ten parameters are an overall 
distance scaling factor i] (relative to an initial arbitrary distance scale) plus the nine first 
and second moments, W and C^. 

Hawley et al. (1986) and Strugnell, Reid, & Murray (1986) were the first ones to 
use a maximum likelihood analysis to determine these 10 parameters. Actually both 
studies use 11 parameter models of the form introduced by Murray (1983), but hold fixed 
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the 11th parameter, the dispersion in rj. As we show in §3.1, current data do not allow 
one to put any useful constraint on this parameter, and we remove it from our primary 
analysis. The dispersion in rj also gives rise to Malmquist (1920) bias. However, in the 
analysis conducted by Hawley et al. (1986) and Strugnell et al. (1986), the inclusion of the 
dispersion in rj does not by itself correct for Malmquist bias. Correcting for the Malmquist 
bias requires additional modeling of the selection criteria of the stars in the sample and of 
the distribution of absolute magnitudes. It is most effectively done on a star by star basis 
(for a discussion see e.g. Smith 1987, also Ratnatunga & Upgren 1997). However when 
the detailed selection criteria are not available and even the magnitude of the scatter is 
uncertain, such procedures are no better than correcting the final result for the "average" 
Malmquist bias (e.g. Jung 1970). 

Ideally, the statistical parallax method should be applied to a group of stars that are: 

1. dynamically homogeneous, i.e. all the stars are drawn from a single velocity 
distribution (not constrained to be Gaussian) regardless of their locations 

2. standard candles, i.e. all have the same absolute magnitude. 

Theoretically, condition (1) can be met by careful selection of stars in the nearby solar 
neighborhood. In practice, to obtain a statistically satisfactory sample often requires that 
stars be collected from a region comparable (~ 3 — 5 smaller) in size to the distance to 
Galactic center. In such a case one does not expect condition (1) to be exactly satisfied. 
However, one can reasonably assume that in this part of the Galaxy the bulk motion and 
velocity ellipsoid change monotonically with distance from Galactic center. Hence, it is 

2 Actually, for small scatter in absolute magnitude, the two approaches are equivalent 
because they differ only by negligible second order terms. 
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important to probe parts of the sky towards and away from the Galactic center. Although 
it is difficult to check whether condition (1) is satisfied, it is possible to investigate some 
of the possible systematic errors analytically and to quantify others through Monte Carlo 
simulations. Condition (2) must be relaxed somewhat, and it suffices to have stars with 
a small scatter around calibrated absolute magnitudes, as we discussed above and as we 
further discuss in §3. 

The RR Lyrae stars of the Galactic spheroid constitute a very good sample for analysis 
by statistical parallax. They are close to satisfying condition (1), being abundant in the 
solar vicinity and distributed reasonably evenly in the sky. Their absolute magnitudes can 
be obtained from the absolute magnitude - metallicity relation (e.g., Carney et al. 1992) 
with a scatter not likely to exceed 0.15 magnitudes (condition 2). Additionally, RR Lyrae 
stars can be observed not only in the Milky Way, but also in neighboring galaxies like the 
LMC. 

The organization of this paper is as follows. In §2 we analyze a simple model of a 
stellar system with an isotropic velocity ellipsoid and we estimate the error in the distance 
scaling parameter rj. We also discuss the role of observational errors and the systematic 
error caused by the mis-estimation of these observational errors. In §3, we apply the results 
of §2 to investigate various possible systematic effects that may bias the statistical parallax 
solution. In §4, we analyze the full maximum likelihood formulation of the problem for 
the limiting case of negligible measurement errors. This case mimics the actual situation 
because the intrinsic velocity dispersion, being much larger than the observational errors, 
dominates the statistical uncertainty. We obtain algebraic expressions for the uncertainty 
in all 10 parameters, rj, Wi, and CV,-. In particular, the error in r\ for the general case has 
the same form as it does for the naive model of §2. That is, as we anticipate in §2, the 
seemingly complicated "black box" of maximum likelihood can be understood in simple 
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physical terms. In §5 we conduct the complete maximum likelihood analysis of the general 
case including observational errors in stellar velocities. We then reanalyze the Layden et al. 
(1996) sample of 162 halo stars found in the solar neighborhood and confirm most of their 
findings. In §6 we conduct Monte Carlo simulations to confirm our analytic results and find 
possible biases, e.g., those induced by the anisotropic spatial distribution of the stars in the 
sample. In §7 we analyze a sample of halo stars constructed by combining the samples of 
Layden et al. (1996) and Beers & Sommer-Larsen (1995). Finally, in §8 we summarize our 
results and discuss their implications for the local distance scale. In particular, we conclude 
that the statistical parallax method is very robust against various systematic effects. 



In this section we conduct the heuristic analysis of the statistical parallax method. It 
is aimed at explaining the physical foundations of this powerful method and gives the basis 
for the analysis of possible systematic errors. In §4 and §5, we develop a more rigorous 
mathematical treatment which is then used in our Monte Carlo simulations described in 
§6. As mentioned in §1, in the general case there are 10 parameters of the fit: the distance 
scaling parameter rj, 3 components of bulk motion and 6 independent components 
of symmetric velocity covariance matrix C^. We initially consider a stellar system with 
an isotropic velocity ellipsoid. The number of independent components of is thereby 
reduced from 6 to 1: CV,- = cr 2 5ij. 

We define the scaling parameter 77 by 



where L true is the actual luminosity of the chosen class of star (e.g., RR Lyrae), and 
L assumed is its assumed luminosity, which can be chosen arbitrarily for purposes of making 



2. Isotropic velocity dispersion 
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the calculation. 

Below we estimate analytically the error associated with the determination of the 
distance scaling parameter 77. When dealing with errors we shall apply the general rule that 
combining errors of statistically independent quantities we add their variances to get the 
total variance, computing the resultant error of two statistically independent measurements 
of the same quantity, we add variances harmonically. 

We decompose the velocity of each star into its radial and tangential components, 
u = (u r , u t ). Now we can conveniently separate two sources of information about r\. The 
first comes from forcing equality between the velocity dispersions as determined from 
radial and transverse velocities ("classical statistical parallax"). The second comes from 
forcing equality in the bulk motion inferred from radial and transverse velocities ( "secular 
parallax"). We assume that these two estimates of 77 are independent in the statistical 
sense. To simplify the analysis we initially assume that there is no bulk motion, which 
implies that (u) = 0, where the brackets, here and afterwards, symbolize averaging unless 
otherwise noted. Thus 



where u t = \u t \. The inferred transverse speed u t is related to the proper motion \i of a star 



"a p ^assumed p j 

7] 7] 

where d trU e and d assume d are the true distance and the distance inferred based on the 
assumed luminosity, respectively. 




(2) 



by: 
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We combine equations (0) and (|3]) to express r\ as 
V 2 = M (4) 



Each component of the velocity, Uj, of a given star may be regarded as being drawn 
from a one-dimensional distribution f(ui). The fractional error in the velocity dispersion 
based on a sample of N stars is 



A(uf) \var(u 2 )} 2 /K-l\* 



(«?> " («?> v iV 

where in the last step we have used the definition of the kurtosis K for (ui) = 0: 



(5) 



* = 7% < 6 > 



By equation 



i 



A(uj) (K-iy 

(u 2 r ) \ N J ■ UJ 

Recall that in equation (Hp designates the sum of the squares of the transverse velocities 
in two perpendicular directions. Thus 

Adding the fractional errors from formulae (0) and @ in quadrature gives the fractional 
error in r\ (see eq. (fD), 



Ar] 
V 



'3{K- 1)Y 2 Gaussian ) / 3 \ 3 
8iV ) ViiV/ 



(9) 



disp V 

where we have used Arj/r) = (l/2)A(r) 2 )/r] 2 , and where in the last step we have evaluated 
the expression for the case of a Gaussian distribution. 
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The bulk motion of the whole sample relative to the Sun gives independent information 
about rj. For simplicity we assumed in the above analysis that the bulk motion was zero, 



where the index i indicates the component and the indices r and t mean "as inferred 
from radial and transverse velocities", respectively. In general, we decompose the bulk 
motion into components in an arbitrary frame of reference which generally results in all 
three components having non-zero values. Here for the sake of simplicity, we pick the 
frame of reference with z-axis aligned with the bulk motion of the sample. We restrict 
attention to the radial velocities of the stars but the results are representative of each of the 
components. One can express the radial velocity as the bulk motion along the line of sight 
with uncertainty equal to the velocity dispersion a: 

u r = W 3<r cos 9 ± a (11) 

where W 3>r represents the bulk motion component inferred from the radial velocity. We 
may construct an estimator of any of the bulk motion components (in this simplified case 
of W 3>r ) by dividing equation (JTTJ) (or its generalized form) by the appropriate angular 
dependence: 



but now we relax this assumption. In analogy to equation (f|) we can write 




Wi 



(10) 



W 3 , r 



a 



(12) 



cos 6 



cos 6 



For the sample of N stars randomly distributed in the sky, the fractional error in estimating 
W 3 , r is 
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where in the last step we have used the fact that (cos 2 #) = 1/3 averaged over all the angles. 
The transverse velocities contain twice as much information as the radial ones, so 



AW 3 ,t = / 3 y a 
W 3tt ~ \2NJ W 3tt 



(14) 



Applying equation (|T0|) and adding the errors from equations (|13|) and ([14]) in quadrature, 
we find that 



At] _ ( 9 \2 a 
~ = V2NJ W 3 



(15) 



A similar analysis can be made to extract information from W\ and W-i- Since these three 
pieces of information about r\ are independent, we add variances harmonically and obtain 



At] 
V 



balk 



(-) 

\2NJ 



W 



K 



K = 



a 



(16) 



where k is the "Mach number", the ratio of bulk motion to dispersion. Combining 
equations (H) and (|l^) yields 



— ' = N'2 

V 



jaussian 



3(K 



-K 

9 



' N~2 [- + -K 2 

3 



(17) 



The typical values for population II stars in the solar vicinity are W ~ 200kms _1 and 
a ~ 100 kms -1 , k ~ 2 so (2/9)k 2 ~ 8/9. Hence, 60% of the information about the distance 
scale comes from velocity dispersions and 40% from the bulk motion. 



2.1. Measurement Errors 

Here we investigate how measurement errors affect the statistical accuracy. We adopt 
the simplified model discussed above with isotropic dispersion a, bulk motion W, typical 
measurement errors a T and <r M , and typical stellar distances D. We define 5 r = a r /a 
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and 5^ = Da^/a. If 77 is determined solely from the bulk motion, then the accuracy is 
(Ar]/r]) 2 = 3A^- 1 k- 2 [(1 + 5 2 r ) + 0.5(1 + On the other hand if 77 is determined from the 
dispersion, (Ary/r?) 2 = 0.5iV- 1 ([(l + 5 2 ) + 0.5(1 + 5 2 )]. Thus, the accuracy of the overall 
measurement is degraded by a fractional amount, 

*<*,/„). *? + S.o.04, (18) 



rv 



(A77/77) 3 6 

where we have adopted <5 r ~ 0.2 and 5^ ~ 0.4 appropriate for RR Lyrae stars. Hence, 
realistic measurement errors have an extremely small effect on the precision of the luminosity 
determination, a result which we confirm numerically in §6. This justifies our approach of 
ignoring measurement errors in our analytic investigation of various other effects conducted 
in §3. 

Most generally, the uncertainty in the statistical parallax determination is set by: 

1. the intrinsic velocity dispersion of the system 

2. the bulk motion of the sample 

3. the observational errors 

4. the size of the sample 

Once the stellar system is given we have no control over items 1 & 2, but we can still vary 
items 3 & 4 to minimize the uncertainty. In the specific case of negligible observational 
errors (meaning small relative to intrinsic velocity dispersion as in the considered case), the 
only remaining factor influencing the uncertainty is the size of the sample. That is, the only 
way to substantially reduce the uncertainty is to increase the sample size and not decrease 
the observational errors. 
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2.2. 



Mis-estimate of the Observational Errors 



As shown in §2.1, observational errors increase only slightly the uncertainty in the 
estimate of rj. Observational errors can also result in a biased estimate of rj if they are 
not properly accounted for. Suppose that the radial velocity errors are ay, but that they 
are taken to be zero in the analysis. (For simplicity, we initially assume that there are 
no proper- motion errors.) If r] is estimated by comparing the radial and proper-motion 
dispersions, the radial dispersion will be overestimated by a factor 1 + (a r /a) 2 while the 
proper-motion dispersion will be properly estimated. Hence rj will be overestimated by a 
fraction (Srj/rj) = [1 + (ay/ 'a) 2 } 1 / 2 — 1 ~ (a r /a) 2 /2. On the other hand, there will be no effect 
on the estimate of r\ based on comparing the radial-velocity and proper-motion bulk motion. 
Hence the total systematic error in the combined determination will be the weighted 
average of the two: (Srj/rj) ~ (a r /a) 2 /[2 + k, 2 /3] ~ 0.3(a r /a) 2 . If the true observational 
error is a r , but the analysis incorrectly treats the error as ay, then (Sr]/r]) ~ 0.35cr 2 /a 2 , 
where 5a 2 = a 2 — a 2 . Finally, we allow for a mis-estimate of the proper- motion errors and 
define 5a 2 = a 2 — a 2 in analogy to the radial-velocity term. We take the typical distance of 
stars in the sample to be D and estimate a net systematic error of 



Since a r /a ~ 0.2 and Da^/a ~ 0.4, the most likely source of a major effect is mis-estimation 
of the proper motion errors, but even this is not likely to be large. For most RR Lyrae stars 
in the sample, the proper-motion error is estimated to be 5 mas yr _1 . Suppose that the 
true value is 20% higher. Then 77 would have been underestimated by only (Srj/rj) ~ 1.7%. 



5i] 



7) 



0.3 



K - d 2 K 

a 2 



(19) 



3. Analytic Investigation of Systematic Effects 
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3.1. Dispersion of the Luminosity 

Throughout our treatment, we have assumed and will assume in §4 and §5 that the 
entire population of stars has exactly the same luminosity. This is not customary. One can 
in principle fit for the dispersion in luminosity (in which case there are 11 fit parameters 
instead of 10) or one can assume a certain dispersion which then affects the values of the 
remaining 10 parameters (e.g. Hawley et al. 1986; Strugnell et al. 1986; Layden et al. 
1996). Here we show that for the RR Lyrae sample, there are several orders of magnitude 
too little information to determine the dispersion from the data as was also suggested 
by previous studies. Moreover, we show that the effect of including the dispersion is one 
order of magnitude smaller than the statistical errors, and that one is therefore justified in 
accounting for this effect separately. 

We assume that the stars have a range of luminosities which we parameterize as a 
dispersion in absolute magnitude M v , o M , but that the analysis is conducted assuming 
that all stars have the same luminosity. We work within the simplified framework 
of §2. The inferred distance to each star in the best-fit solution will then deviate 
from the true distance by 0(om\ This will increase both the variance in the inferred 
transverse speeds and the mean. Averaged over all directions, the fractional increase of 
the dispersion is 2(lnl0/5) 2 [l + K, 2 /6]a 2 M . Thus, if r) were determined by matching the 
dispersions of the radial velocities and proper motions, it would be underestimated by 
a fraction (Sr]/r]) = — (lnl0/5) 2 [l + K 2 /Q]a 2 M ~ —QAa 2 M . On the other hand, if i] were 
determined by matching the bulk motions, it would be underestimated by a fraction 
(Srj/rj) = — 0.5(ln lO/S) 2 ^ ~ — O.lcr 2 ^. Thus, by comparing the two estimates, one could 
in principle determine <jm- For the RR Lyrae sample, however, the precision of each 
method is only ~ 7%, so that the precision of the difference is ~ 10%. This means that 
the data set would have to be increased ~ 400 fold in order to detect plausible values of 
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&m < 0.15 mag. If the dispersion is simply ignored and the two methods are averaged (as 
they automatically are in the maximum likelihood approach), then r\ will be underestimated 

by (5r]/r]) = -(lnl0/5) 2 [12 + 3k 2 ]ctm/[ 12 + 2^ 2 ] 0.25a 2 M . This is an extremely small 

correction for plausible values of o~m- The underestimate in rj in turn induces changes of the 
bulk motion and dispersions, which we estimate numerically in §6. 

3.2. Malmquist bias 

Malmquist bias is another correction to t] which scales cx a 2 ^, but it is of opposite 
sign. As we discuss in the introduction, the Malmquist bias correction is not automatically 
taken into account by simply considering the 11-parameter model, with the 11th parameter 
being the dispersion in 77. The assumption intrinsic to most such models (and also to 
our 10-parameter model) is that the stars in the sample are representative of the whole 
population. This assumption is not correct for samples that have been selected based on 
stellar magnitudes. Malmquist bias must therefore be put by handQ. Malmquist (1920) 
showed that the mean distance of a magnitude-limited sample drawn from a homogeneous 
distribution with a Gaussian scatter in absolute magnitude will be higher than that of one 
without scatter by (Stj/tj) ~ 3 (In 10/5) 2 cx 2 M ~ (In 10/5)1. 38^ ~ +0.64a 2 M . 

Malmquist bias is a particular form of "selection bias" which arises whenever a sample 
deviates systematically from the underlying population because of the way that it was 
selected. The magnitude of the selection bias depends, in general, on the form of the 
selection. The RR Lyrae sample is "local" in the sense that it contains stars that are less 
than 2 or 3 kpc from the Sun, but does not contain substantially more distant stars (even 

3 See e.g. Smith (1987) for alternative approach, but see also the discussion of its limited 
application in §1. 
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though such stars are plentiful). Some procedure (whose details are not known to us) had 
to be applied to select this sample. In principle, one could distinguish between nearby and 
distant stars by several methods, including apparent magnitude, trigonometric parallax, 
and proper motion. However, while we do not know exactly how the sample was selected, 
we do know that of these types of information only apparent magnitudes were available. 
Therefore, each star can have been included only because it satisfied some magnitude limit, 
even if this magnitude limit is different from star to star. The correction calculated by 
Malmquist (1920) is independent of the magnitude limit, so the same correction applies to 
all sub-samples of fixed magnitude limit, even samples with only one object. Therefore, 
it applies to the sample as a whole. For an inhomogeneous underlying population, the 
numerical pre-factor in 3(ln lQ/h) 2 a\ 4 may in general differ from 3. However, we have 
confirmed numerically that for the distribution of RR Lyrae stars from which the Layden 
et al. (1996) sample is drawn, the pre-factor is in fact very close to 3. 



3.3. Velocity-Position Correlations 

Statistical parallax makes the implicit assumption that the stars seen in all directions 
have the same velocity distribution. Since the RR Lyrae stars in the sample have typical 
distances ~ 2kpc, which is a significant fraction of i?o ~ 8kpc, it is possible that this 
assumption is not valid. Suppose, for example, that the velocity dispersion in the z 
direction falls with distance from the Galactic plane. If one attempted to fix the luminosity 
by matching the z-dispersions of the radial velocities and the proper motions, one would 
underestimate its value. This is because the radial-velocity measurements would be made 
primarily on stars far from the plane, while the proper-motion measurements would be 
made primarily on stars near the plane. The high dispersion of the latter would be 
mistakenly interpreted as indicating that the stars were closer and hence less luminous 
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than they actually are. While this systematic bias would be diluted by unbiased dispersion 
measurements in the other two directions and by unbiased measurements of the bulk 
motion, it would still affect the final result. (We focus attention on the z direction because 
that is the only axis for which there is a physical plane of symmetry. Any correlations along 
the other two axes would tend to have opposite signs in the positive and negative directions 
and hence would cancel to lowest order.) 

To make a quantitative estimate of the size of this effect, we consider a simple model 
where the stars are observed over a uniform sphere of radius D, and have a z velocity 
dispersion, o 2 z (z) = a(l — f\z\/D). Consider an ensemble of stars at a distance r and angle 
9 from the north Galactic pole. These will enter into the radial-velocity and proper-motion 
estimates of the z dispersion with statistical weights cos 2 9 and sin 2 9, respectively. If the 
correct luminosity of the stars were adopted, the respective dispersion estimates would be 



Thus, one would be led to make an underestimate (5r)/r)) = —(9/64)/. Since the effect 
is diluted by measurements of the dispersion in the other two directions and of the bulk 
motion, the net bias is (Srj/rj) = — (3/64)//[l + k 2 /Q] ~ —0.03/. For plausible values of 
f<D/R ~ 0.25, this would be a small but perhaps not completely negligible correction. 

We have therefore reanalyzed the sample including an 11th parameter, the velocity- 
dispersion gradient in the z direction. The best fit scale length for this gradient is 
7kpc, but since the reduction in x 2 is om y ^X 2 — 1-15, the gradient is not statistically 
significant. In any event, even if the gradient is real it introduces a systematic error of only 




(20) 



and 




(21) 
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(Srj/ri) = -0.005. 



4. Analytic Predictions 

In this section we use the maximum likelihood method to analyze the errors intrinsic 
to the statistical parallax method. We again neglect velocity measurement errors. We 
argued in §2.1 that for the halo stars of the Galaxy the accuracy of the determination of 
the model parameters scales with the size of the sample of stars and is barely influenced by 
observational errors. Having selected the group of stars (e.g., RR Lyrae stars), we choose an 
orthogonal frame of reference and resolve each star's velocity into components designated 
as Ui. We assume that the velocity distribution is a three-dimensional Gaussian, but we 
do not assume isotropy. At first Gaussianity of the assumed velocity distribution seems to 
be a rather restrictive condition, especially in view of the fact that the actual distribution 
is highly non-Gaussian. However, we argue at the end of this section (§4.1), and confirm 
numerically in §6 that no bias is introduced by assuming that the distribution is Gaussian. 

The probability of finding a star with three velocity components in the ranges 
(ui, Ui + dui) is given by: 



L(ui] r), dj, Wi)d 3 u = — r exp --Y\ (m - Wi)(C - Wj 

assumed (2ll\C\)^ L 2 l 'J 3 ' 



d 3 u, (22) 



where W is the bulk motion relative to the Sun, C is the velocity covariance matrix, and 
\C\ is its determinant. The volume element in three-dimensional true velocity phase space 
is d 3 u, whereas d 3 u\ assumed is the corresponding element in "assumed" velocity phase space. 
Since two of the "assumed" components are proportional to i] -1 , d 3 u = rf ' d 3 u\ assumed . 
Hence the probability density of finding a star with the observed velocity components in 
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the ranges (ui, Ui + dui) is given by: 



L{ui\ 77, Cij, Wi 



(2n\C\Y 



exp 



(23) 



The logarithm of the total probability of finding the system in the observed state is: 



ln£ = ]T^ =i lnL fc 



(24) 



where is the probability ( p3|) associated with the k th star. The curvature matrix is given 
by 



B, 



a/3 



d 2 ln£ 

dad/3 



t-^n d 2 In L k 



(25) 



where a and /3 range over the 10 parameters of the model. The covariances among the 
parameters are then given by 



cov 



a/3 



(26) 



We designate the inverse of the covariance velocity matrix as 



Q = C 



-1 



(27) 



which is a special case of the general definition (|52|). We express 



Ui = U r Ti + T)U, 



(28) 



in terms of its radial and transverse components. Note that t{ is the projection of u t , 
defined in §2, on z-th direction. We also define Vi to be the random part of the velocity: 



Vi = Ui- Wi. 



(29) 
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We evaluate the curvature matrix in two steps. First we evaluate the contribution of 
each star to equation (p4]). Next, we sum over all stars under the assumption that they are 
isotropically distributed on the sky. We find 

d 2 In L 1 

-^QlmQnk + {Qv) n Qlm{Qv) k + (k <-+ I) + (m Tl) (30) 



dC mn dCki "2 
d 2 InL 



drjdCki 

^ InL 

3W s dC kl 



(Qt) k (Qv)i + (k~l) (31) 



Qsk(Qv)i + {k^l) (32) 



<9 2 lnL 2 

-^r = ^ + mt) (33) 
92lUL -(Qt)s (34) 



dW s drj 

^ InL 

dWidW a 



Qsi, (35) 



where the symbol (k «-> /) means "add the terms with k and / exchanged, but only if k ^ V , 
and where we have employed Dirac notation, i.e. 

(X|(9|y)=V XiOM. (36) 

Averaging over positions leads to: 
1 d 2 In £ 1 

QimQnk + (k <-> Z) + (m <-> n) (37) 



NdC mn dC kl 2 
1 <9 2 ln£ 2 



Qki + {k^l) (38) 
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1 d 2 In £ 

= (39) 



NW a dC M 
1 «9 2 ln£ 



jy 9r/ 2 ?7 2 3?7 2 
1 d 2 \nC 2 



^ + (40) 



(41) 



N dW s dr] 3r] 

1 <9 2 ln£ „ fAn . 

~ — ^ — = Qsi (42) 

NdWidW a K 1 

The inverse of the matrix B, given by equations ( j37"D - (f42"D , describes the variances 
and covariances of the parameter estimates (see eq. fl26D). Note that no observational errors 
have been included in this treatment and that the uncertainties in the model parameters 
scale with the square root of the size of the stellar sample. We now adopt the frame of 
reference for which the axes are aligned with the axes of the velocity ellipsoid i.e. where CV, 
is diagonal. We then find that the errors in rj, Wi, and are given by 



var{r\) a 

7] 2 = N 

var(Wi) 1^,4^ 2 



(43) 

QuWfa) (44) 



C u N V 9 

Var(Gjj) ~ * r (l + la) (45) 

(46) 



C] 3 .W 9 



var(Cij) 1 
CuCjj N 



where 

,3 9 

Here k is the generalization of the Mach number to an anisotropic distribution. 



a=(^ + lK 2 ) k 2 = (W\Q\W). (47) 
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We have derived fully analytic estimates of the errors in all parameters. Note that the 
result given in equation (f43[) can be written as 



That is, the error in the distance estimate (eqs. ( |4"8| ) and (f43|)) for the general case has the 



same form as the error for an isotropic distribution derived using a few basic assumptions 
(eq. (0)). 



4.1. Non-Gaussian Velocities 

As we show in §6, the actual distribution of RR Lyrae stars is very far from Gaussian. 
On the other hand, the formalism presented in §4 explicitly assumes that the velocities 
are Gaussian. That is, the likelihood given in equations (^) and fl2"3|) is the probability 
distribution for a 3-dimensional Gaussian. At first sight this appears to be a very serious 
problem because the true form of the distribution is only crudely determined from the 
data and there is no a priori argument by which one knows even how to parameterize the 
distribution. We address this problem in two ways. First we argue in this section that 
one does not introduce a bias by using a Gaussian likelihood function regardless of the 
form of the actual velocity distribution. (The use of a Gaussian function does cause one to 
incorrectly estimate the uncertainties of the luminosity measurement, but by a calculable 
and - as it turns out - small amount.) Second we confirm this result by Monte Carlo 
simulations in §6. 

Why does the assumption of a Gaussian distribution not introduce a bias? Early 
versions of the statistical parallax methods were explicitly based on the first and second 
moments of the velocity distribution and as such were intrinsically insensitive to the form 
of the stellar distribution. In the maximum likelihood method we should ideally choose 
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the velocity distribution function that matches the actual one. For the case where there is 
not enough information about the underlying velocity distribution, it is desirable to pick 
a Gaussian. This is because the Gaussian likelihood has the special characteristic that it 
works in essence by determining the means and dispersions of the velocities (in each of 
three dimensions) separately from the radial velocity and proper-motion measurements and 
then forcing these to be equal by fixing the luminosity. The means and dispersions (and 
also covariances) are then reported as nine additional parameters of the fit. The method 
therefore effectively reproduces the naive procedure outlined in §2, which was based on the 
first and second moments of the distribution. That is why it also reproduces the results 
derived using that procedure. It is straightforward to show that if one uses maximum 
likelihood to fit a non-Gaussian distribution to a Gaussian function parameterized by mean 
W and variance a 2 , then the resulting values of W and a 2 will be unbiased estimators of the 
mean and variance of the (non-Gaussian) distribution^. Thus, adjusting the luminosity to 
maximize the Gaussian likelihood of a non-Gaussian distribution still amounts to equating 
the means and variances of the (non- Gaussian) distribution. Since the determinations of 
these means and variances are unbiased, so is the estimate of the luminosity. 

If the underlying distribution is non-Gaussian, the Gaussian maximum likelihood 
procedure will return estimates of the errors that differ from the true errors. Consider 
for example the simple isotropic model with dispersion a, kurtosis K, bulk motion W, 
and negligible observational errors. The maximum likelihood estimate of the error in 
rj is (see eq. ([43])) (Arj/i])^ 2 = (2/9)iV[6 + k 2 ] while the true error is (see eq. (|17|)) 
(Ai]/r])~ 2 = (2/9)N[(12/[K — 1]) + k 2 ]. For a more general velocity distribution that is the 
product of distributions in the n, 8, and z directions with dispersions (a n , ag, a z ) and kurtoses 

4 To be more precise a 2 will be an unbiased estimator of (N — 1)/Nx the variance, but 
this slight difference has no practical impact on the discussion here. 



-24- 



(K n , Kg, K z ) (and assuming that the bulk motion is in the 6 direction), equation ( f!3"D 
yields (Ar]/r])~ 2 = (2/9)7V[6 + (W/ag) 2 ]. Generalizing from equation fllTD , we estimate the 
true errors as (Ar//r/)- 2 = (2/9)N[(A/(K 7T - 1) + A/(K e - 1) + A/(K Z - 1) + (W/a e ) 2 }. 
(We confirm this estimate numerically and mention some practical complications in §6.) 
For an arbitrary velocity distribution, the true and maximum likelihood-estimated errors 
could in principle differ substantially. However, for the actual RR Lyrae population, 
4/ (K w — 1) + 4/ (Kg — 1) + 4/ (K z — 1) ~ 7, close to the Gaussian value of 6. This implies 
that the maximum likelihood estimated errors are nearly equal to the true errors. 



5. Complete analysis 

In this section we obtain the formulae needed to analyze real data and to carry out 
Monte Carlo simulations. Our analysis here is similar to the one conducted in §4, but we 
account for a few additional effects. 

The probability density of finding a star with velocity components in the ranges 
(uj, Ui + duj) is now given by: 



rf r 1 

L(ui\ 77, Cij, Wi) = - exp --^2^( Si - Wi )(M ) ij (s j - wj) 



(2n\M\y- 



(49) 



where Wi is the bulk motion (defined more precisely below), and Si is the stellar velocity 
expressed in its local Galactic coordinate frame (the velocity that would be measured by an 
observer located at the star's position and at rest with respect to the Galactic center) 



S i = J2 j=1 R H ( u i ~ v ®i) ' ( 50 ) 
with v = (-9,232,7) kms" 1 . 
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The matrix M is the full velocity covariance matrix, defined as 



(51) 



Here a r is the observational error in the radial velocity, o t is observational error in each of 
the inferred transverse velocity components, and P^- is the projection operator on the plane 



velocity. The velocity v is the sum of the velocity of the Sun's LSR vlsr and the peculiar 
motion of the Sun relative to its LSR v Qj lsr- The completely new element not present in 
formula (|23"D (nor in any previous statistical parallax analyses of which we are aware) is 
the rotation operator i?^ that accounts for the fact that the star's velocity is drawn from 
a certain distribution (e.g., three-dimensional Gaussian) in the Galactic frame of reference 
centered at the star which is in general a very different distribution when expressed in terms 
of the Galactic frame centered at the Sun. The matrix elements of the rotation operator are 
functions of the star's position and the distance from the Sun to Galactic center. We shall 
assume that the distance to Galactic center is ~ 8kpc. The bulk motion, Wi, is relative 
to the local Galactic frame. That is, in the Galactic frame at the Sun's position, the bulk 
motion is J2j(R~ 1 )ijWj. As in the no-error case, the logarithm of the total probability of 
finding the stellar system in the observed state is the product of probability functions of 
the form (|49|) for all the stars: ln£ = J2k=i l n -^fc- 

We generalize equation (p7|) to 



of the sky. In equation ( pUD we explicitly subtract the velocity of the Sun v 



from the star's 



Q = M 



-i 



(52) 



We maintain the decomposition ttj = u r ri + r/ti, and we introduce the short notation 




(53) 
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As previously, we use Dirac notation ( |36D to designate scalar products and the symbol 
(k <-> I) means "add the terms with k and I exchanged, but only if k ^ /". 

The first derivatives of InL are: 

<91n L 2 

drf = ~v + ~ W ^ X \ P ^ + ^ a t tT (.Q P ) ( 54 ) 

dlnL . . 

= -Xi (55) 



dWi 
dlnL 1 



dC kl 2 



(Qki ~ x k xi) + (k<r+l) (56) 



Second differentiation leads to: 



^ = ^ + (y\Q\y)-a 2 (x\P\x)-A W 2 (y\QP\x) + a 2 tr(QP) 



d 2 InL 

dwidwj 



-2r] 2 aftr(QPQP) + 4ri 2 a?(x\PQP\x) (57) 
Q lj (58) 



d 2 In L 1 

QlmXnXk - -QlmQnk + (/ <"> fc) + ijn <-> 7l) (59) 



dCkldC mn " 2 

<9 2 InL 



drjdwi 
d 2 InL 
dr]dC H 
d 2 InL 



dwjdC, 



hi 



-{Qy)i + 2 w 2 (QPx) t (60) 
-:r fe (Qy), - W 2 (QPQ) M + 2r l a 2 (QPx) k x l + (/<->*) (61) 
Q«2* + (/ «-> fc) (62) 



In deriving these results we have made use of the identities 
<91n|M| / ,0Af\ dM" 1 ,9M„ , 

We apply equations (53) - (^2j) and Newton's method to find the maximum likelihood 
solution for the Layden et al. (1996) sample of halo RR Lyrae's. Table 1 gives the basic 
results and illustrates the process of correcting them for detected biases. Our "assumed" 
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distance scale is the one derived by Layden et al. (1996). As a result, if our determination 
were in perfect agreement with Layden et al. (1996), the scaling parameter r\ should be 
exactly equal to 1. The four rows give the maximum likelihood values of the ten parameters 
characterizing the sample under four different sets of assumptions. None of the rows 
corresponds exactly to the results of Layden et al. (1996), but the first row with rj increased 
by ~ 0.015 does so approximately. The first row is not corrected for rotation or biases. 
The values presented in the second row are obtained using a non-unit rotation operator 
(see discussion following eq. (|5"T|)). Including rotation increases r] by about 0.4%, only 

i 

slightly influences the bulk motion, but makes the velocity distribution less triaxial (C 2 2 2 

i i 
closer to C 3 2 3 ) and more elongated (higher C^). The third row gives our best estimate for 

the values of the parameters after correcting for the biases discussed in §6. We also take 

into account that maximum likelihood underestimates the velocity variances by a factor 

N/(N — 1) as mentioned in §4.1. Finally the fourth row takes into account the scatter in 

the absolute magnitudes of RR Lyrae stars. The direct effects are estimated numerically 

by introducing artificial scatter in the Monte Carlo simulations and are in good agreement 

with the estimate given in §3.1. The indirect effects of Malmquist bias are as described 

analytically in §3.2. 

The stars that enter our analysis are those defined as "Halo-3" population from Table 
3 in Layden et al. (1996). Also, we use Layden's best estimate of the proper motion errors 
(6.5 mas yr -1 ) for stars taken from the catalog compiled by Wan, Mao & Ji (1980) rather 
than the value (5masyr _1 ) artificially adopted by Layden et al. (1996) to match their own 
sample (although this makes almost no difference). We assume that RR Lyrae stars follow 
the absolute magnitude-metallicity relation (My = const + 0.15[Fe/H]) of Carney, Storm & 
Jones (1992). However, the results are only sensitive to the value of the absolute magnitude 
at the mean metallicity of the sample, ([Fe/H]) = —1.61, and not to the slope of the relation. 
We checked that the solutions for different slopes are statistically indistinguishable from one 
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another. By taking account of the metallicity, we restrict the possible scatter in RR Lyrae 
absolute magnitudes to the intrinsic scatter at fixed metallicity. To make the corrections 
discussed in §3, we adopt for the Layden et al. sample of field Galactic RR Lyrae stars 
om = 0.15. To make this estimate, we first inspected color-magnitude diagrams of several 
globular clusters and found that for these very homogeneous populations the dispersion 
is typically o~m — 0.08. We take this as a lower limit. The only available inhomogeneous 
sample at approximately fixed distance is the 28 RRab Lyrae stars of Hazen & Nemec 
(1992), which are 4° east of the LMC center and for which a M ~ 0.17. We take this as 
an upper limit because some of the scatter may be due to the dispersion in distance. For 
example, if the LMC RR Lyrae population has an r -3 5 profile and the flattening c/a = 0.6 
(similar to the Galactic population), we find a dispersion due to distance CM,dist = 0.09, 
implying an intrinsic dispersion of <7m ~ 0.14. In any event, those prefering other values 
of o"m should note that it is straightforward to find the corrected values of all parameters 
simply by scaling the difference between rows 3 and 4. The correction we apply in Table 1 is 

^ « 0.25a 2 M - 0Ma 2 M = -0.39a 2 M . (64) 

where the first term, estimated numerically, is due to the scatter in the absolute magnitude 
of RR Lyrae stars (§3.1) and the second term is due to Malmquist bias (§3.2). 

To make connection to previous studies of Galactic structure, note that it is customary 
to use the following notation: 

ill 
(Tu = C±\, o~y = C 2 2 2 , cr w = C33. (65) 

Notice, however, that by using the non-unit rotational operator R we are actually measuring 
the underlying velocity distribution in the local Galactic frames of the stars in the sample 
under the assumption of Galactic axisymmetry. We therefore use (ir, 9, z) rather than 
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(U, V, W). The rotation operator R is the appropriate first order correction to the rectilinear 
solution (e.g., Layden et al. 1996) regardless of whether the velocity distribution in the 
Galaxy is exactly axisymmetric. Our best estimate of the RR Lyrae absolute magnitude at 
the mean metallicity of the sample ([Fe/H]) = -1.61 is My = 0.75 ± 0.13. The velocity 
ellipsoid is (a w , a e , a z ) = (171 ± 10, 98 ± 8, 96 ± 8) kms -1 and the RR Lyrae population is 
moving in 9 direction at —211 ± 12 kms -1 relative to the Sun, where the errors are adopted 
from YYTA case in the Table 3. 



6. Monte Carlo simulations 

Here we present the results of Monte Carlo simulations aimed at checking and 
fine-tuning our analytic results. There are two main classes of simulations. One class 
keeps all the star positions from the Layden et al. (1996) sample unchanged giving results 
closely related to those obtained using the real sample. For the second class, the stars are 
placed randomly over the celestial sphere. The position-related biases can be determined 
by comparing the two classes of simulations. 

In each simulation we construct 4000 mock samples of 162 halo stars (162 is the 
size of the Layden et al. 1996 sample). For each sample we generate a set of 162 space 
velocities drawn from a distribution with specified means, dispersions and kurtoses in the 
three principal directions. For each star we transform its velocity components from the 
star's Galactic frame of reference to the Sun's frame and, in some cases, add Gaussian 
measurement errors in accordance with the values given by Layden et al. (1996). 

In the second step we find the most probable parameters describing each of the samples. 
To analyze our mock samples we use exactly the same maximum likelihood procedure that 
was used to obtain the results for the real stars. 
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We perform several tests to check for the various systematic effects. The results of 

these investigations are summarized in the Tables 2 and 3. Table 2 gives the biases in 

i 

rj, Wi, and Cj found in each simulation. The first column says whether the observational 
errors were included or were set equal to zero. The second column tells whether the 
velocities were rotationally adjusted (operator R in equation fl50l) ). The third column 
gives the information on whether the positions of the stars in the sample were the true 
ones from Layden et al. (1996) or were chosen randomly. The fourth column lists the 
assumed values of kurtoses in all three directions^]. The fifth column assigns the name to 
each case. The name contains the most important information about the case. For example 
YYTH means that the observational errors were included in the analysis, that velocities 
were rotationally adjusted, that the true positions of the stars were considered, and that 
the kurtoses in all directions were higher than Gaussian (here equal to 4). Generally, "T" 
stands for true, "R" for random, "G" for Gaussian, "H" for high, "L" for low and "A" for 
actual. The next seven columns give the results for the most interesting parameters of the 
fit based on 4000 realizations. In all cases, the biases in the off-diagonal elements of the 
velocity covariance matrix (normalized in the same way as in the Table 1) are smaller than 
0.01 and we therefore do not display them here. The input values of the parameters for 
the underlying distribution are given below the descriptions of the column content. The 
values in parentheses in the first row give the approximate errors in the determinations of 
the biases based on the NNRG case. The values in the table are the biases B detected 
for a given case defined as: B = (obtained value) — (template value). The exceptions to 
this rule are the biases of the dispersion which are computed according to the scheme: 
B = (obtained value) — [(N — l)/N] 5 (template value). This comes from the fact that 

5 Note that for YYTA case, the input kurtoses are the ones that produce (in the mean) 
the same output values as those obtained for the actual sample. 
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maximum likelihood returns variances that are (N — l)/N times the values of the true ones 
(see §4.1). Note that with this definition of bias, one must subtract B from the maximum 
likelihood solution to get the corrected value. 

The first five rows test the correctness of our implementation of the maximum likelihood 
method. They contain the cases with Gaussian velocity distributions (K w = K e = K z = 3) 
for which we expect the biases to be small. Indeed, the deviations from the input values in 
the cases with the random positions of the stars (NNRG and YYRG) are extremely small, 
that is, the obtained solutions are statistically indistinguishable from the input. For the 
cases with the true positions of the stars (NNTG, YNTG and YYTG) the deviations remain 
small but now they are statistically significant and should be treated as biases rather than 
statistical fluctuations. 

The next three rows, below the blank line, give the results for the samples with 
velocities drawn from high- or low-kurtosis, non-Gaussian distributions. One may suspect 
that our maximum likelihood procedure derived assuming a Gaussian velocity ellipsoid will 
not be very accurate in such cases. To the contrary, the determination of the parameters is 
extremely robust and insensitive to the underlying velocity distribution. The case YYTA is 
constructed to reproduce the input kurtoses of the real sample of the Layden et al. (1996) 
stars. A glance at Figure 1 shows that the distributions of the stars in the sample are highly 
non-Gaussian. The measured kurtoses of RR Lyrae stars in three principal directions are: 
(2.18,3.14,3.93). The right panel of Figure 1 shows how close the distribution in the radial 
direction is to a box distribution. Analyzing our simulations we find that the most probable 
values of the kurtoses of the true underlying distribution of RR Lyrae stars in the solar 
neighborhood are (K n , Kg, K z ) = (2.04,3.22,4.28). Fortunately, as we see from the lower 
part of the Table 2, the parameter determination is almost completely independent of the 
form of velocity distribution. 
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Table 3 contains the errors in the parameter determinations as obtained from the 
scatter of the Monte Carlo realizations. For comparison the first row gives the errors 
predicted analytically by equations (|43"D - (f£6l) . For the simulation with velocities with 
no measurement errors and unit matrix operator R (NNRG) our analytic estimates of the 
errors should agree with the scatter of the simulated results. The agreement is striking and 
all the discrepancies are within statistical uncertainties. Additionally, accounting for the 
true positions of the stars (NNTG), observational errors (YNTG) and non-unit rotation 
operator (YYTG) affects the errors only slightly. Furthermore, the same is true for the 
high-kurtosis (YYTH) or low-kurtosis (YYTL) cases. The errors in the non-Gaussian cases 
cannot be predicted exactly from the formulae given in §4.1 because the distribution is not 
isotropic. Nevertheless, the predicted trend is clearly confirmed. We conclude that the 
analytic estimates of the uncertainties constitute excellent approximations to the realistic 
case. Both observational errors and higher kurtosis slightly increase the uncertainties. 

In Table 4, we compare the kurtoses measured in the three principal directions to those 
of the underlying velocity distributions. First note that unless the underlying kurtosis 
is very low (as in the NYRL case), the measured kurtosis is underestimated due to the 
finite size of the sample (e.g., NNRG and NYRH cases). Additionally, measurement errors 
tend to Gaussianize the distribution, and as a result reduce kurtoses in high-kurtosis cases 
(e.g., YYRH vs. NYRH), increase kurtoses in low-kurtosis cases (e.g., YYRL vs. NYRL), 
and have little effect on a Gaussian distribution (e.g., YYRG vs. NNRG). The measured 
kurtoses are also biased due to the non-isotropic positions of the stars in the sample on the 
sky (e.g., NNTG vs. NNRG or YYTH vs. YYRH), but this effect is relatively less important 
for non-Gaussian distributions. Our best estimate for the underlying kurtoses of the halo 
RR Lyrae stars, based on the simulation results presented in Table 4, is (2.04, 3.22, 4.28) 
(YYTA). 
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7. Kinematics of Metal-Poor Halo Stars 

As we emphasized in §2.1, the factor that fundamentally limits the precision is the 
size of the sample. This limitation is severe because to expand the sample requires first 
identifying the new RR Lyrae stars and then measuring their proper motions. Some 
additional stars could undoubtedly be found in the South with distances similar to those 
characteristic of the sample. In most cases, however, new proper motion studies would be 
required. Moreover, the increase in the sample size would be modest. To really increase 
the size of the sample substantially would require finding stars at greater distances. In the 
North, at least, proper motions would be available from the Lick Proper Motion study. 
However, the size of the errors for these stars would begin to approach those of the proper 
motions themselves, increasing the possibility of systematic errors. Thus, this route would 
appear to require substantial additional work, lasting perhaps several decades. In addition, 
stars that are found too far from the Sun may have different kinematics which would 
introduce additional systematic errors. 

There is, however, a poor man's route to increased statistics which has the side-benefit 
of yielding new kinematic information about the stellar halo. If one repeats the analysis of 
§2, but with different numbers of stars with radial velocity (N r ) and proper motion (N t ) 
measurements, one finds 




in place of equation (|T7|). If one measured only the radial velocities for a substantial new 
population of RR Lyrae stars, one could hope to drive down the error by up to a factor 3 1 / 2 
without obtaining any new proper motions. In fact, it is not even necessary that they be RR 
Lyrae stars. The only requirement is that the radial-velocity stars and the proper-motion 
stars have the same kinematics. This is trivially satisfied if all the stars are RR Lyrae, but 
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it can also be attained through careful selection of other stars. 

For this purpose, we turn to the catalog of 1936 non-kinematically selected metal-poor 
([Fe/H]< —0.6) stars of Beers & Sommer-Larsen (BSL, 1995). We select stars both from 
the BSL sample and Layden et al. sample with [Fe/H]< —1.5. This should insure that the 
great majority are from the stellar halo, not the thick disk. (If the fraction of the thick disk 
stars differs between the two samples, this will introduce systematic errors, as we quantify 
below.) Next, we exclude BSL stars classified as "variable" to avoid overlap with the RR 
Lyrae sample and also to avoid stars with possibly poor radial- velocity determinations. 
Finally, we eliminate stars with estimated distances > 3 kpc. The precise distances are not 
important for the analysis, but we do wish to eliminate stars that are outside the volume 
from which the RR Lyrae stars are drawn since they may partake of different kinematics. 
These criteria yield a BSL sub-sample of 724 stars and a Layden et al. sub-sample of 
106 RR Lyrae (of which 103 are from the sample of 162 stars analyzed elsewhere in this 
paper). The BSL stars are overwhelmingly turn-off, giant-branch, upper- main- sequence, 
and horizontal-branch stars (in that order) and so have masses (or progenitor masses) very 
similar to those of RR Lyrae stars. Hence, they should be drawn from very nearly the same 
underlying population, regardless of the process by which the stellar halo formed. 

The new combined non-kinematic sample consists of 830 stars total: 724 stars come 
from the BSL sample and 106 RR Lyrae stars from Layden et al. In Table 5 we compare 
the parameter and error determinations for the combined non-kinematic and Layden et al. 
samples. The results are corrected for magnitude dispersion, Malmquist bias, and other 
effects described in §5 and §6. The Malmquist bias correction for the non-kinematic sample 
is the same as the one applied in Table 1 of §5. The magnitude dispersion correction was 
determined based on numerical simulations. Consequently, we compare final results: the 
best values and uncertainties for the Layden et al. (1996) sample are repeated from the final 
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rows of Tables 1 and 3, and the corresponding data are presented for the non-kinematic 
sample. The two results appear quite consistent, but since they are based partly on the 
same data, a careful statistical analysis is required. 

The first point to note is that Table 5 gives only the statistical error, i.e., the error 
that arises primarily from the finite size of the sample and secondarily from measurement 
errors. For the non-kinematic sample, there is an additional error arising from the possible 
differences in the fraction of thick-disk stars in the two sub-samples. Let these fractions be 
fi and fii respectively. Let f — (fi + /2V2, and let A = (f\ — /2V/- Since the proper 
motions are determined entirely from the Layden et al. sub-sample and the radial velocities 
are determined almost entirely from the BSL sample, the thick disk stars will generate an 
additional error 5rj/r] ~ gfA where (1 — g) ~ 0.5 is the ratio of the typical thick-disk to 
typical halo speeds relative to the Sun. We estimate the thick disk fraction / from the 
fact that the non- kinematic sample has a rotation speed of 35kms _1 (see Table 5). Since 
halo stars have approximately zero rotation (BSL), this net rotation must be produced by 
thick-disk contamination. Since the rotation speed of the thick disk is of order 170kms _1 
(Casertano, Ratnatunga, & Bahcall 1990), the contamination level is / ~ 35/170 ~ 0.2. 
The mean distance from the plane is very similar for the two sub-samples, 1210 pc for the 
Layden et al. stars and 926 pc for the BSL stars. We therefore estimate g ~ ±0.2, and hence 
Srj/rj ~ ±0.02. Since this error is uncorrelated with the statistical error, we add the two in 
quadrature to obtain our final estimate for the non-kinematic sample r\ = 0.953 ± 0.054. 

Next, to either compare or combine the two samples, we must evaluate the correlation 
coefficient of the two determinations, 7 = c^cn^)" 1 / 2 - Here is the covariance matrix 
of the two measurements, with en = (0.057) 2 (from Table 5) and c 2 2 = (0.054) 2 (derived 
in the previous paragraph). In the appendix, we find 7 = 0.46. This is very close to 0.5, 
the value one would naively guess because the radial-velocity measurements are almost 



completely independent while the proper motions are almost completely dependent. 

The 1 a expected difference between the two measurements is therefore 
(en + c 22 - 2ci 2 ) 1/2 = 0.059. That is, the two measurements, 771 = 0.980 ± 0.057 
corresponding to My = 0.75 ±0.13 at [Fe/H]= —1.61, and r] 2 = 0.953 ±0.054 corresponding 
to My = 0.79 ± 0.12 at [Fe/H]= —1.79 differ by less than la. It is therefore appropriate to 
combine them. We apply standard linear theory (e.g., Boutreux & Gould 1996 eqs. 2.5 and 
2.8) to find a best estimate and error 



corresponding to My = 0.77 ± 0.011 at [Fe/H]= -1.71. 

Finally, we note that Table 5 contains the first complete solution of the velocity ellipsoid 
of halo stars as determined from a non-kinematically selected sample. One important 
(though not unexpected) result is that the stellar halo is not moving relative to the LSR 
in either the radial or vertical directions. From Table 5, the estimates would appear to be 
3 ± 8 km s _1 away from the Galactic center and ± 5 kms" 1 toward the north Galactic pole. 
In fact, each of these values (and errors) should be augmented by a factor (1 — f)^ 1 ~ 1.25 
because the contaminating thick disk stars are known to be bound to the disk and so have 
zero mean motion in both directions. Hence the best estimates for the halo bulk motion in 
the two directions are 4 ± 10 km s" 1 (radial) and ± 6 km s" 1 (vertical). The rotation of the 
halo (35 ± 8 kms" 1 ) is larger than in the solution for the 162 RR Lyrae alone, but this is to 
be expected because that sample was selected in part kinematically, that is by eliminating 
stars with the most prograde orbits (Layden et al. 1996). 



7/ 



m - 



(Vi - 772) (en - C12) 
Cll + c 22 - 2Ci2 



0.965 




8. Discussion and Conclusions 



We have investigated many potential sources of systematic error in the RR Lyrae 
absolute magnitude calibration by statistical parallax using a combination of analytic 
and Monte Carlo techniques. We find that all corrections to previous results are small 
and, in particular, that the highly non-Gaussian RR Lyrae velocity distribution does 
not bias the determination at all even though (or rather, because) the method explicitly 
assumes a Gaussian distribution. We find that the mean RR Lyrae absolute magnitude is 
My = 0.75 ± 0.13 at the mean metallicity of the sample ([Fe/H]) = —1.61 compared to 
My = 0.71 ± 0.12 obtained by Layden et al. (1996). The largest source of difference comes 
from including Malmquist bias which makes our estimate 0.03 mag fainter (for our adopted 
scatter a M = 0.15). Most of the rest of the difference (0.01 mag) comes from the other 
corrections for scatter in the absolute magnitude (0.03 mag brighter for Layden et al. versus 
0.01 mag brighter for us). There are several other smaller effects that reduce the difference 
by 0.01 mag. 

We also analyze a semi-independent non-kinematically selected sample of stars 
with metallicities [Fe/H]< —1.5 taken from Layden et al. and BSL and find similarly 
My = 0.79 ± 0.12 at ([Fe/H]) = —1.79. Additionally, this analysis yields measurements of 
the radial bulk motion (4 ± lOkrns -1 ) and vertical bulk motion (0 ± 6kms _1 ) of the halo 
relative to the LSR. 

Our principal result is therefore that the RR Lyrae absolute magnitude calibration by 
statistical parallax is extremely robust and that the statistical error (0.13 mag) should be 
taken at face value. 

If one assumes that type ab RR Lyrae stars in the LMC are (apart from different mean 
metallicity) similar to those in the Layden et al. (1996) sample, then the distance modulus 
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to the LMC is 



/xlmc = (V)oxmc - {0.75 + 0.15([Fe/H] + 1.61)}, 



(68) 



where (V)o,lmc is the dereddened mean apparent magnitude of RRab's at the distance of 
the center of mass of the LMC. Ideally, the way to determine this quantity is to measure 
V for a large sample of RRab's in an annulus around the LMC bar. Stars well away from 
the bar should be little affected by internal LMC extinction and the mean foreground 
extinction is reasonably well understood. Moreover, by taking the average of annulus, one 
would assure that the mean distance of the sample is equal the center-of-mass distance, 
regardless of the geometry of the RR Lyrae distribution. Such photometry should soon 
be available from ongoing microlensing surveys. However, the three RR Lyrae samples 
currently available are both smaller and less ideal. One is Walker's (1992) RR Lyrae sample 
drawn from 6 LMC clusters (excluding one foreground cluster) for which 



The error is from the scatter and hence accounts for measurement errors and uncertainty in 
the mean distance of the clusters relative to the center of mass. We estimate an additional 
error in the mean extinction of 0.03. Applying equation (|6"S|), we find /zlmc — 18.28 ± 0.14. 
However, two lines of evidence suggest that cluster RR Lyrae stars may be systematically 
brighter than those in the field. First, Sweigart (1997) has argued that abundance anomalies 
seen in cluster giants but not in halo field giants could be due to rotational mixing in the 
former. Such mixing would dredge up helium in horizontal branch stars, making them up 
to a few tenths of a magnitude brighter. Second, new main-sequence-fitting distances to 
Galactic globulars (Reid 1997; Gratton et al. 1997) would, if confirmed, imply that cluster 
RR Lyrae stars are several tenths of a magnitude brighter than the value for field stars 
reported here. Hence, it seems prudent to restrict the comparison to LMC field RR Lyrae 
stars. 



(^)o,lmc = 18.98 ±0.03 



«[Fe/H]> 



1.9) 



(clusters). 
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Two photometric studies of LMC field RR Lyrae stars have been conducted, both 
in B. Graham (1977) found (B) = 19.61 ± 0.02 for 60 field RRab stars near NGC 1783, 
about 4° northwest of the bar. The error reflects only the dispersion in the measurements 
and Graham (1977) notes that calibration errors are possible. There is no consensus 
on the reddening toward this field. We adopt the estimate of Alvarado et al. (1995), 
E(B — V) = 0.08 ±0.02. From the period-amplitude diagram of these stars, Hazen & Nemec 
(1992) estimate a mean metallicity ([Fe/H]) = —1.3. To estimate the mean color, we use 
the relation (B - V) = 0.658 + 0.710 log P + 0.097[Fe/H] from Caputo & De Santis (1992) 
which, from their Figure 5b, has an uncertainty of ±0.015 for each individual star. The 
mean period of the Graham (1977) stars is (P) = 0.564 days, implying (B — V) = 0.36. 
Using R v = A V /E(B - V) = 3.1, we then find 

(V) ,N1763 field = 18.92 ±0.08, 

where no account has been made for a possible calibration error. 

Hazen & Nemec (1992) found (B) = 19.61 ± 0.03 for 28 field RRab stars near NGC 
2210, about 4° east of the bar. They believe that their photometry is well calibrated. They 
adopt E(B - V) = 0.08 ± 0.01 from Caldwell & Coulson (1985), and find ([Fe/H]) = -1.8 
from the period-amplitude diagram. The mean period is (P) = 0.576 days, which implies 
(B - V) = 0.31. Hence 

(VOo,N2210 field = 18.98 ±0.05, 

Combining the apparent magnitude and metallicity measurements with the absolute 
magnitude of RR Lyrae stars yields distance-modulus estimates /iNi763 field = 18.18 ± 0.15 
and /iN22iofieid = 18.26 ± 14 which seem quite consistent. However, unless the distribution 
of RR Lyrae is spherical, these fields will not be at the same distance as the LMC center of 
mass. Assuming that (like the Galactic RR Lyrae stars) they have a radial profile oc r~ 3 - 5 
and an axis ratio c/a = 0.6, and taking the inclination of the LMC to be i = 27° (Bessel, 
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Freeman & Wood 1986), we find offsets of +0.05 mag and —0.05 mag for the two fields, 
respectively. This leads to two estimates for the LMC center of mass /iLMc 3field — 18.13±0.15 
and /iLMC° ficld = 18.31 ± 0.14. Since the errors of these determinations are correlated, 
the difference is A/x = 0.18 ± 0.09, which is uncomfortably large. While both of these 
determinations are much smaller than the traditional value /iLMC — 18.5 and are smaller yet 
compared to several new determinations (Reid 1997; Gratton et al. 1997; Feast & Catchpole 
1997), there remains the possibility that the underlying measurements of (V)o are affected 
by systematic errors. A robust estimate of this quantity should therefore await results from 
the microlensing surveys. 
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A. Determination of the Correlation Coefficient 

We now return to the framework of §2 in order to evaluate the correlation coefficient 
7 = Ci2(cnC22)~ 1 / 2 - The error for a single determination of 77 has the form given by 
equation (|66|) . For two partially dependent determinations of rj, this generalizes to 

% = ^(6 + K 2 )- 1 (2cI i + 4) (Al) 

where c^- is the 2x2 covariance matrix of the errors in rj, 

( 1 N r ,12 



V 



N r ,l N r<1 N r , 2 
N r , 12 1 



(A2) 



and similarly for c*. Here N r i = 162 is the number of radial-velocity stars in the Layden 
et al. sample, N r> 2 — 830 is the number in the non-kinematic sample, and N r ^ 2 = 103 
is the number of overlap stars. Similarly, N t> i = 162, N t> 2 = 106, and N t ^2 = 103. It is 
immediately clear that the diagonal elements of equation (|A1|) reduce to equation (|66|). 
The off-diagonal terms take account of the covariance between the two radial-velocity-based 
velocity ellipsoids and the covariance between the two tangential-velocity-based velocity 
ellipsoids. However, equation (|A1|) takes account only of the statistical errors. We now 



take account of the additional systematic error caused by differering levels of thick disk 
contamination. Recall that this error affects only the non-kinematic determination, so it is 
included by adding a matrix to equation (|A1|) 



c ij = |(6 + ^Y\2c% + c%) + (gfAf I ° ° I , (A3) 



where gfA = 0.02 was evaluated in §7. Hence, 7 = 012(011022) 1 ' 2 = 0.46. 
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Fig. 1. — Left panel shows the velocity distribution functions for three Galactic directions. 
Right panel shows that the distribution in radial direction is very boxy. 
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Table 3. Errors for the cases considered in our Monte Carlo simulations 



I I I 

Case name A77 Awj Aw 2 Aw 3 AC& AC 2 \ AC 3 2 3 AC 12 AC 13 AQ 





0. 


.052 


13. 


.20 


10 


.67 


7 


.55 


10. 


,98 


6.53 


6 


.27 


0, 


.079 


0. 


.079 


0, 


.079 


NNRG 


0. 


.052 


12, 


.95 


10 


.64 


7. 


.61 


10. 


99 


6.47 


6 


.19 


0, 


.079 


0. 


.079 


0, 


.080 


NNTG 


0. 


.053 


12. 


.87 


11 


.26 


7. 


.63 


11. 


62 


6.69 


5. 


.79 


0, 


.078 


0. 


.080 


0, 


.080 


YNTG 


0. 


.057 


13. 


.26 


11 


.84 


8. 


.08 


12. 


.14 


7.24 


6 


.29 


0. 


,086 


0. 


.086 


0, 


091 


YYTG 


0. 


.056 


13. 


.25 


11 


.73 


8. 


.08 


12. 


09 


7.21 


6 


.28 


0. 


.085 


0. 


.086 


0, 


090 


YYRG 





,054 


13. 


.31 


10 


.97 


8. 


,16 


11. 


.45 


7.01 


6 


.81 


0. 


.085 


0. 


.087 


0, 


.091 


YYTH 





,059 


13. 


.28 


11 


.96 


8. 


,12 


13. 


.98 


8.26 


7. 


.39 


0, 


.086 


0. 


.086 


0, 


.091 


YYTL 





.055 


13. 


.60 


11 


.62 


8. 


.12 


10. 


20 


6.20 


5 


.12 


0. 


.088 


0. 


.086 


0, 


.089 


YYTA 





.057 


13. 


.59 


11 


.83 


8. 


.14 


10. 


37 


7.50 


7. 


.64 


0. 


.086 


0. 


,084 


0, 


.090 



Note. - As previously, the name of the case uniquely characterizes statistical properties 
of velocity input. The first row, just below column description, gives the errors predicted 
analytically in §4. 



Table 4. Comparison of measured kurtoses to those of underlying velocity distribution. 
The effects of the finite size of the sample, observational errors and non-isotropic positions 

of the stars in the sample are presented 



Case name Underlying kurtoses Measured kurtoses 



NNRG 


(3,3,3) 


(2.98, 


2.99, 2.98) 


YYRG 


(3,3,3) 


(2.99, 


3.05, 3.05) 


NNTG 


(3,3,3) 


(2.99, 


2.97, 3.01) 


YYTG 


(3,3,3) 


(3.01, 


3.02, 3.05) 


NYRL 


(2,2,2) 


(2.02, 


2.03, 2.03) 


YYRL 


(2,2,2) 


(2.14, 


2.34, 2.36) 


YYTL 


(2,2,2) 


(2.16, 


2.36, 2.31) 


NYRH 


(4,4,4) 


(3.91, 


3.91, 3.91) 


YYRH 


(4,4,4) 


(3.82, 


3.73, 3.72) 


YYTH 


(4,4,4) 


(3.83, 


3.66, 3.77) 


YYTA 


(2.04,3.22,4.28) 


(2.18, 


3.14, 3.93) 



Data from Layden et al. (1996) 




-400 -200 200 400 -400 -200 200 400 
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TABLE 1 

Value for the ten parameters for Layden et al. (1996) sample of halo RR Lyrae stars 















i 


i 








Description 


V 




w 2 


w 3 






C33 




Cl3 


C23 


No rotation, raw output 


0.9808 


-0.14 


23.70 


-4.51 


166.45 


102.41 


95.82 


-0.1019 


0.0471 


-0.1420 


Rotation, raw output 


0.9844 


0.53 


23.26 


-4.62 


169.38 


98.85 


96.04 


-0.0975 


0.0472 


-0.1302 


Rotation, "numerical" biases corrected 


0.9888 


0.97 


21.02 


-5.49 


171.15 


98.64 


96.18 


-0.0954 


0.0505 


-0.1416 


Rotation, correction for My scatter 


0.9801 


0.95 


20.55 


-5.55 


171.27 


98.11 


96.10 


-0.0952 


0.0502 


-0.1422 



Note. — The bulk motion and velocity dispersions are given in kms 1 and the off-diagonal coefficients of the velocity 
correlation matrix are shown in normalized dimensionless form: CV, = Cij(CuCjj)~* . 



TABLE 2 

The biases in 8 cases considered in our Monte Carlo simulations 



Observational Rotation Positions Kurtoses Case r\ W\ w 2 w 3 C 2 2 2 C 3 2 2 

errors adjusted of stars (K x ,K y ,K z ) name 1.0 0.0 22.0 -5.0 168 100 96 



No 


No 


Random 


(3,3,3) 


NNRG 


No 


No 


True 


(3,3,3) 


NNTG 


Yes 


No 


True 


(3,3,3) 


YNTG 


Yes 


Yes 


True 


(3,3,3) 


YYTG 


Yes 


Yes 


Random 


(3,3,3) 


YYRG 


Yes 


Yes 


True 


(4,4,4) 


YYTH 


Yes 


Yes 


True 


(2,2,2) 


YYTL 


Yes 


Yes 


True 


(2.04,3.22,4.28) 


YYTA 



(0.0008) 


(0.20) 


(o. 


,17) 


(0, 


,12) 


(0 


.17) 


(0 


,10) 


(0.10) 


0.0008 


-0.25 





,13 


0, 


17 


-0 


,16 


-0, 


02 


-0.07 


-0.0040 


-0.52 


1. 


,96 


0, 


.78 


-1 


,12 


0, 


46 


0.08 


-0.0042 


-0.54 


2, 


,18 


0, 


,86 


-1 


,23 


0, 


.54 


0.14 


-0.0042 


-0.45 


2, 


,20 


0, 


.87 


-1 


,25 


0, 


55 


0.17 


0.0012 


-0.27 





,12 


0, 


.15 


-0 


,11 


-0, 


13 


-0.15 


-0.0040 


-0.47 


2 


,19 


0, 


.88 


-1 


,21 





,47 


0.21 


-0.0028 


-0.44 


1. 


,83 


0, 


63 


-1 


.03 


0. 


,54 


0.33 


-0.0044 


-0.44 


2 


,24 


0, 


,87 


-1 


,24 


0. 


,51 


0.16 



Note. — The input values of the parameters for the underlying distribution of stars are given below the descriptions of the 
column content. The values in parentheses in the first row give the approximate errors in the determinations of the biases based 
on the NNRG case. 



TABLE 5 

Comparison between two samples 

























Description 


V 












^33 




Cl3 




Layden et al. sample 


0.980 


0.95 


20.55 


-5.55 


171.27 


98.11 


96.10 


-0.095 


0.050 


-0.142 




(0.057) 


(13.59) 


(11.83) 


(8.14) 


(10.37) 


(7.50) 


(7.64) 


(0.086) 


(0.084) 


(0.090) 


BSL and Layden ([Fe/H]< -1.5) 


0.953 


2.85 


34.91 


0.37 


161.65 


108.59 


94.43 


0.014 


-0.087 


-0.014 




(0.050) 


(8.30) 


(8.49) 


(4.65) 


(6.51) 


(7.60) 


(4.95) 


(0.083) 


(0.066) 


(0.085) 



Note. — The bulk motion and velocity dispersions are given in km s 1 and the off-diagonal coefficients of the velocity correlation 
matrix are shown in normalized dimensionless form: CV, = Cij(CuCjj)~* . Quantities in parantheses are the uncertainties in 
parameter determinations. 



